library(Seurat)
library(tidyseurat)
library(data.table)
library(harmony)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')

# read in matrix data ----------
cd45p <- fread('esophagealCancer/data/GSE160269_CD45pos_UMIs.txt.gz')

cd45p[1:5, 1:5]

# check if rownames are unique
cd45p <- column_to_rownames(cd45p, 'V1') |>
  as.matrix() |>
  as('dgCMatrix')

sobj <- cd45p |>
  CreateSeuratObject(names.delim = '-',
                           min.cells = 3,
                           min.features = 200)

# modify metadata ----------
sobj$orig.ident |> table()

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |>
  VlnPlot('mito.ratio', pt.size = 0)

# 111028 > 109983
sobj <- sobj |>
  filter(mito.ratio < 10)

# remove normal tissue
sobj <- sobj |>
  filter(str_detect(orig.ident, 'T'))

ITlist <- read_tsv('esophagealCancer/data/ITlist.txt',col_names = 'id')

sobj <- sobj |>
  mutate(genotype = case_when(
    orig.ident %in% c('P8T', 'P63T', 'P130T') ~ 'TT',
    orig.ident %in% ITlist$id ~ 'IT',
    .default = 'II'
  ))

sobj |>
  VlnPlot('mito.ratio', pt.size = 0, group.by = 'genotype')

sobj <- sobj |>
  mutate(barcode = .cell) |>
  separate(barcode, into = c('sample', 'immune', 'cell'))

sobj |> count(immune)

flt_meta %>%
  mutate(barcode2 = case_when(
    str_detect(foo, '126|127|128|130') ~ barcode,
    TRUE ~ str_remove(barcode, 'T')
  )) -> flt_meta

# read from excel ------------
meta_from_paper <- readxl::excel_sheets('data/source_ESM.xlsx') %>%
  set_names() %>%
  map(read_excel,
      path = 'data/source_ESM.xlsx',
      skip = 1)

pub_cell_meta <- meta_from_paper[["Metadata for fibroblasts"]] %>%
  select(c('cell', 'celltype', 'tissue', 'sample')) %>%
  mutate(barcode2 = str_replace_all(cell, '\\.', '-'))

pub_cell_meta %>%
  separate(barcode, into = c('foo', 'bar','tag'), remove = 0) -> pub_cell_meta

left_join(flt_meta, pub_cell_meta, by = 'barcode2') -> bind_meta

bind_meta %>%
  column_to_rownames('barcode') -> sobj@meta.data

sobj <- subset(sobj, genotype != 'Normal')

sobj <- SetIdent(sobj, value = 'genotype')

# pre-process --------
sobj <- sobj |>
  quick_process_seurat()

# rds checkpoint =====
sobj |> write_rds('esophagealCancer/zhang2021escc.rds')

sobj <- read_rds('esophagealCancer/zhang2021escc.rds')

DEenrichRPlot(sobj, 
              ident.1 = 'NAF2',
              enrich.database = 'GO_Biological_Process_2021',
              max.genes = 3000,
              return.gene.list = 1) -> delist

DimPlot(sobj, group.by = 'genotype')

# find markers for clusters
harmonySobj %>%
  FindAllMarkers(only.pos = TRUE,
                 min.pct = 0.1) %>%
  group_by(cluster) %>%
  top_n(n = -10, wt = avg_log2FC) ->
  topMarkers

write_csv(topMarkers, 'results/walker_ESCC_TopMarker.csv')

write_rds(harmonySobj, 'data/walker_ESCC.rds')

# walker dataset read in
load('data/walker_OAC.DGE.RData')

OAC.DGE[1:5, 1:5]

CreateSeuratObject(OAC.DGE,
                   names.delim = '_',
                   min.cells = 10,
                   min.features = 3) -> sobj

sobj@assays$RNA[1:5,1:5]

sample <- read_tsv('data/walker_sample.txt')

sample_idents <- sample$sample

sobj <- subset(sobj, idents = sample_idents)

VlnPlot(sobj, 'nCount_RNA', log = TRUE)

sobj <- subset(sobj, nCount_RNA > 200)

metadata <- sobj@meta.data

metadata$genotype <- 'II'

metadata$genotype[which(metadata$orig.ident %in% c('OAC1411T', 'OAC174T'))] <- 'IT'

metadata$genotype[which(metadata$orig.ident == 'OAC132T')] <- 'TT'

ggplot(metadata)+
  geom_bar(aes(x = genotype, fill = genotype))

write_csv(metadata, 'results/walker_meta.csv')

sobj@meta.data <- metadata

SetIdent(sobj, value = 'genotype') -> sobj

load('data/walker_Cell.Metrics.RData')
write.csv(Metrics, 'results/walker_pub_meta.csv')